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Abstract — ^We demonstrate that it is possible to compute wave 
function normalization constants for a class of Schrodinger type 
equations by an algorithm which scales linearly (in the number of 
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1. Introduction 

IN a recent paper[l] it was demonstrated that it is possible to 
solve some eigenvalue problems of the Schrodinger type to 
almost arbitrary high precision. The cases presented explicitly 
were (i) the ground state eigenenergy of the anharmonic 
oscillator, 

-V'" + a;^V' = eV', (1) 

which was found to an accuracy of more than 1 000 000 
decimal digits, (ii) the eigenstate number 50 000 of the same 
equation where the eigenenergy was found to an accuracy of 
more than 50 000 decimal digits, and (iii) the lowest even and 
odd parity states of the double-well potential, 

-s2v>:^ + (a;2-l)2^±=e±V±, (2) 

where the two eigenenergies were found to an accuracy of 
more than 30000 decimal digits for the case s = 1/50 000. 

In the latter case the two eigenenergies are degenerate to 
almost 29 000 decimals, with the difference directly com- 
putable by the WKB method. The 10* order WKB expansion 
given in [2] provide an accuracy of about 48 decimals for the 
difference, all of which agrees with the difference between the 
two numerical calculated eigenenergies. 

Also the eigenenergies of the highly excited states of 
equation (1) can be calculated by the WKB method. The 12* 
order WKB expansion given in [3] provide an accuracy of 
about 67 decimals, all of which agrees with the numerical 
result. 

It is certainly difficult to find physical systems where one 
need to know eigenenergies to tens of thousands of decimals 
or more. However, if one need to compute the wavefunction 
very accurately to very high values of x (to f.i. evaluate 
matrix elements) the value depends extremely sensitively on 
the eigenvalue parameter. Thus, in our opinion, algorithms for 
very-high-precision evaluation of eigenvalues may be useful in 
combination with routines for very-high-precision evaluation 
of matrix elements and normalization integrals. 

In this paper we demonstrate that the latter can be achieved 
very simply, with the number of eigenfunction evaluations 
growing only linearly with the desired precision P in dec- 
imal digits. Each eigenfunction evaluation will in turn re- 
quire a number of high-precision multiplications which grows 
linearly with P, and each such multipUcation require a 
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CPU time which scales asymptotically between P^'^ and 
P log P log log P (depending on which high-precision mul- 
tiplication algorithm is used). Thus, the total time to evaluate 
one normalization integral to P decimals precision can be 
expected to increase somewhat faster than P^. However, since 
the eigenfunction evaluations are independent they can be run 
in parallel. 

The rest of this paper is organized as follows. In section II 
we make some general remarks on numerical integration 
rules. We base our analysis on the Euler-Maclaurin suimnation 
formula and the Poisson resummation formula. Our conclusion 
is that the trapezoidal rule is not only the simplest one but also 
the best one. For integrals over finite intervals some there are 
endpoint corrections which should be considered separately; 
these corrections vanish for our normalization integrals. In 
section III we consider some simple example cases similar 
to our wavefunction normalization integrals. For an infinite 
(in principle) integration range and a fixed number M of 
integration steps one must strive for a balance between the 
error due to using a finite stepsize h, and the error e(.Tmax) 
due to covering only a finite integration range. These errors can 
be estimated by respectively analysing the Fourier transform 
f{p) of the integrand f{x) as \p\ — > oo (by the method of 
steepest descent), and the asymptotic behaviour of f{x) as 
I a; I 00. It appears that both analyses can be extended 
to wavefunction normalization integrals by use of the WKB 
approximation. This extension is done in section IV. 

II. Remarks on integration formulae 

There is no scarcity of numerical integration formulae in 
the literature [4]. Standard choices are rules for /{/} = 

fl^dx f{x) which reproduces integrals of polynomials below 
a certain order exactly, 

I{.f}^li.fo + fi)h^T{h), 

I{.f}^-A.h+ifi+f2)h = S{h), 

9 (3) 
Hf}^l {/o + 3/i + 3/2 + h)h = Ss/s{h), 

^if}-J^ (^/o + + 12/2 + 32/3 + 7/4) h = B{h), 
45 

with h = {b — a)/M and /„ = /(o + mh) when there 
is M -I- 1 terms in the integration rule. These are known 
respectively as the trapezoidal, Simpson's, Simpson's |, and 
Boole's rule. They are automatically exact for polynomials 
which are antisymmetric about the midpoint x = {b — a)/2. 
The M independent coefficients are chosen to give exact 
results for all symmetric (about x) polynomials of order below 
2M. However, as M increases the weight coefficients develop 
in a suspicious way. To integrate functions over large intervals 
to very high precision it seems dubious to extend the procedure 
above. 



A. Extended Simpson 's rule 



rules 



An alternative way to handle integration over large intervals 
is to divide them into many smaller ones, and apply one of the 
rules above to each of the subintervals. We denote this by a 
bar over the rule, T{h) — > T{h) etc. A rather common choice 
is to extend Simpson's rule, leading to the formula 



J a 



dxf{x) « S{h) =\{h + 4/1 + 2/2 + • • • 

+ 2/M-2 + 4/M-l+/A/)/i, 



(4) 



which looks curious to any member of an equal society. What 
is wrong with half of the points? Are they of the wrong 
gender? Which mysterious force of numerical error analysis 
leads to this spontaneous breakdown of translation invariance? 
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Fig. 1 Comparison of accuracy of the extended trapezoidal and 
Simpson's rules for some integrands with very smooth boundary 
behaviour. 



In our opinion only T{h) is a logically sensible procedure. 
We claim that this is also the best procedure. To support and 
exemplify this claim we have evaluated the integrals below: 



Ax (1 - a;2) 



2\" 



Joo = / da; 



1 — tanh 



2a; - 1 
1 - (2a;- 1)^ 



(5) 
(6) 



by the T{h) and S{h) = ^T{h) - \T[2h) rules. The relative 
errors of the two methods as function of discretization length h 
are shown in figure 1. The result is as expected that S behave 
as well as T with twice the discretization length h. 



B. Euler-Maclaurin summation formula 

The examples in figure 1 are untypical in the sense that the 
integrands are very smooth at the endpoints (which is typical 
for the normalization integrals we are primarily interested in). 
One may use the Euler-Maclaurin summation formula to see 
which quantity is actually computed by the various integration 
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where A^''^ = [,P'''>{b) - f''\a)] h^+^, and the ellipsis in 
the first line denote "non-perturbative" terms depending on 
how f{x) varies in the full integration range. With S{h) — 
T{2h) one finds the corresponding quantity for the 
impson's rule: 
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I.e., the A*^i^ -coefficient has been eliminated at the cost of 
increasing all the remaining coefficients'. 

A natural way to eliminate the A *^i) -coefficient would be 
to calculate /' explicitly at the endpoints, or use a numerical 
approximations like 
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+ -f'^-'>ix)5^ + . 



(3) 
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for sufficiently small 5. One should take 5 somewhat smaller 
than h. The modified trapezoidal rule described in [5] is 
generally worse than the extended Simpson's rule because of 
an inaccurate approximation of /'. However, there are few 
reasons to use the same discretization length for computing 
derivatives at the endpoints as for computing the bulk contri- 
bution to the integral. To avoid introduction of significant new 
errors it is sufficient to choose 5 < h/2 in the first formula. 
By choosing S = h/\^20 in the last two one also eliminates 
the A^'^-'-term in equation (7). 



C. Poisson resummation formula 

There exist expressions for the "non-perturbative" terms in 
equation (7), but we find the Poisson resummation formula 
more clarifying: Assume f{x) is an entire, integrable function, 
with Fourier transform 

/(p)= r dxf{x)e^\ (9) 

J —00 

Then one version of the Poisson resummation formula reads 

(10) 
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This can be used to estimate the accuracy of the numerical 
integration formula 



dxf{x) 



00 



^ hfimh)-J2f{'^-^). (11) 

m= — 00 fc^O 



When f{x) is an entire function its Fourier transform f{p) 
will vanish faster that any inverse power of p as p — >^ cx). This 
means that the numerical approximation to the (infinite range) 

'Note that the error terms given for the "extended" rules in reference [4] 
are in disagreement with our results. 



integral of an entire function will converge very fast towards The saddle point equation becomes 

the exact value as /i — > 0. jj i \ _ n 2n-i • _ n 

Integrals of f{x) over a finite range [a,b] can be viewed ^ {^stP) = ~'^n'Xg +vp = \j, (19) 

as integrals of g{x) = 9{b - x) e{x - a) f{x) over an ^jth 2n - 1 solutions 
infinite range. But g(x) will usually be discontinuous, with 

a fourier transform g(p) which vanishes only algebraically as ^ _ ^m{4k+i)/{4n-2) / ^/(^"-i) _q ^ 2n — 2 

p 00. The "perturbative" endpoint corrections of the Euler- * \2nJ ' > 

Maclaurin formula can be used to account for these algebraic (20) 

terms in a systematic way. But there is no point in making The two most relevant saddle points^ occur for fc = and 

endpoint corrections beyond the error in the bulk contribution, n - 1, at which the real part of the exponent (l>{x,p) is (when 

the latter being of magnitude /(^) + f{-^). P ^ '^'^ 



III. Example INTEGRALS 



Re^(..,p) = -(2n-l)sin(^) (^) 



2n/(2n-l) 



In this section we will analyze the behaviour of some simple = anh ^"/(^" (21) 

cases which are similar to typical ground state normalization t,,,. , ^. ...^ ., 

integrals leadmg error due to a finite stepsize h is of magnitude 

gRe0(a;a,p) _ g-a„h 2n/(2n i) ^-^p^j. ^ prefactor of less 

2 importance). On the other hand, including only the M first 

A.- e ^ terms of the sum leads to an error of magnitude e~K^+^)''l " . 

Consider first the ground state of the harmonic oscillator, For a given M we should choose h to balance these errors, 

/(a;) = e-^', in which case i-^- (i_i/2„) 



/>) = ^/^e-f l\ (12) 

Hence we have 



6„ (M + l)^^'-'/""\ (22) 

(2n-l)/4n^ 



/oo 
-oo 



with &„ = (^)'/'" [(2n-l)sin(j^) 
leads to an error of magnitude 



This 



da;e--' = ^ /ie-('"'*)'-V4;^e--'/'''--- . (13) e(M) = e'^-^^+i), (23) 



In practise we must approximate the infinite sum by a finite 
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(2n-l)sin(5;f^) 
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I.e., to evaluate 



one, the integral numerically to P decimals accuracy one must 

choose 



M 



dxe-^' =/i + 2/i5]]e-(™'')' M+1 > i^^Pw (0.12 + 0.467 n)P. (24) 

+ 2/ie^*^^+^^^''^ + • • • — \/47re^'^^/''^ _ . . . where the numerical approximation is very good for n > 2. 



If we make h, too small the first correction term (on the second q g^-i^^-a'^ f 

line above) becomes loo large. If we make h too large the sec- ' / 2 2\2 

ond correction term becomes too large. I.e., the optimal choice Finally consider f{x) — ~°- ) . In this case 
of h occurs approximately when e~*^^+^^ = e~'^ 1^ , 



7(a) = / 

J — C 



This leads to an expected error of order j 2-ir (1/4) as a ^ 

e(M) = e"''^''^"'"^^ (15) ~^ \ V^/a as a ^ 00. 

I.e., to evaluate this integral numerically to P decimals accu- ^he saddle point approximation to the Fourier transform 
racy one must choose 



(25) 



/OO P'OO 
da;e-(^'-"')'+'f^ = / daje-^^^-f), (26) 
-00 J —00 



M + 1 > ^°^Pw0.73P. (16) 

TT leads to the saddle point equation 



B. e-^ 

In more generality consider f{x) = e~^ with n a positive We introduce 
integer. In this case the integral is 

Kn= j dxe-"'" = . (17) 



(j)'{xs,p) = -xl + a^Xs+i^ =0. (27) 



3y 



2n I ' which leads to a quadratic equation for y^. 



We use the saddle point method to estimate its Fourier ^3 _ .j. (29) 

transform 4 27y^ 

/oo /.oo 
jj2.g-a'^"+ipa: = / da;e'^^^'^^. (18) ^The deformation of the integration path so that it passes through these 
^ 7-00 saddle points is an interesting exercise[6] 



It is convenient to rewrite p so that (p/8) ~ (a/\/3) sinhSr;. 
Then the solutions for equation (29) becomes 
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,±(r,+i7r/6) 



(30) 



I.e.. 



±(7;+i7r/6) 



3y 

which in both cases of ± leads to the solution 

2a 



V3 



cosh {rj + in/G) . 



(31) 
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Fig. 2 Predicted (lines) and obtained (points) precision as function of 
M, for the three example integrals discussed in the text. In order from 
top to bottom: / Axe~'^'^ ~° ' > / Axe'"^ , J da^e"'' . The number 
of function evaluations in each case is A/ + 1. As can be seen, the 
obtained precision agrees well with the theoretical estimate. 

There are two more ways to take the cube root. One of them 
leads to an equally relevant saddle point, 



2a 

71 



cosh {rj + i57r/6) 



(32) 



while the last one is irrelevant. The real part of the exponent 
at the relevant saddle points is 

Re0(a;s,7?) = — -a* sinh^ rj cosh 277. (33) 

This provides an error estimate in parametric form: If we 
choose a finite stepsize 



(34) 
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Aa^ sinhSr/' 
the corresponding error will be of magnitude 

£{h) ~ g— f a" sinh^ )) cosh 2r; 



(35) 



The error caused by summing over only a finite range 
of x-values, < x^i„ < x < a;,Tiax, should be cho- 
sen to be of the same magnitude as e{h). I.e., with s = 



(4/3) sinh^ 77 cosh 2ri, 



= a(l + ^/i)^/^ 

if s > 1, 



•^min 



(1 — ^/s)^^'^ Otherwise. 



(36) 
(37) 



I.e., we must use M = (xmax — x^inj/h evaluation steps in 
the numerical integration. 

IV. Wavefunction normalization integrals 

Our investigation of the example integrals gives us con- 
fidence that we can obtain a fairly good a priori estimate 
of the obtainable precision e{M) at a given number M of 
discretization steps, at least asymptotically for large M. For 
this we need to (i) estimate the behaviour of the Fourier 

transform, 4''^{p), of the integrand at large p = 211 jh (to find 
the obtainable accuracy e(K) at a given stepsize A), and (ii) 
estimate how the integrand decays away from its maxima (to 
find the required x-range of summation for the same accuracy). 

Both quantities can be obtained to reasonable accuracy by 
use of the WKB approximation. 

A. WKB estimates of wavefunctions for a:^" potentials 

For large x an estimate of solutions to the eigenvalue 
problems 

-ij}" + (a;2" - 1/^ = 0, 
can be written in the form^ 



(38) 
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where Xq" — E, and 



C = exp 
— exp 
= exp 



n + l 



dt- 
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1 n~l 
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2(n + l) 

TT / TT \ 
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2 \2nJ 
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1 ^(«+l)/2n 



N 



(39) 



with B{x,y) = T{x)r{y)/r{x + y) the Beta function. We 
have evaluated the WKB integral using partial integration, and 
let X — ?> 00 in the remainder. Further, in the last equality 
we have assumed E = Ej^ to be eigenvalue number N , and 
evaluated it by the WKB approximation. 

We use this approximation to estimate the Fourier transform. 



EE / dxiP{xfeP'', 



(40) 



by the saddle point method. The saddle point equation can be 
written 



= E-{p/2)\ 



(41) 



givmg 

ip'^{p) 



exp(^Re {ipx,)) 

e.^{~^,s^{^)p[{p/2f-Ef'-]. 

(42) 



'We ignore the algebraic prefactor (x^" — E)' 



-1/4 



Thus, to compute the normahzation integral to P decimals 
precision one must choose p ~ 2Ti/h so that the right hand 
side of (42) becomes equal to 10^ (or smaller). The value 
of h is best found numerically. But note that one at least must 
have {p/2f > E, or 

For large x the normalization integrand behaves like 

%l){xf « exp ^-2 dt\/t^" - 

« exp (^-^xVx^'^-E^ . (43) 

Thus, to compute the normalization integral to P decimals 
precision one should integrate to a value x — x^ax for which 
the right hand side of (43) becomes less than 10^^. The value 
of Xjnax is best found numerically. But note that one at least 
must have x^^ > E, or 

pl/2n 

This provides a lower limit on the number of integration steps. 



M + 1 = > i £'(»+l)/2« 

h TT 



/_jdWi 



(44) 



Here we have assumed E = Ej^j to be eigenvalue number N, 
and used the WKB approximation to estimate its value. Since 
N is the number of nodes in the eigenfunction, the inequality 
(44) quite reasonably says that the number of evaluation points 
must be larger than the number of oscillations. 



B. Ground state wavefunction of x potential 

Specializing the results of the previous section to the ground 
state of the a;^ -potential, and approximating Eq by zero, one 
finds from equation (42) that 



and from equation (43) that 

i^ixf 



(45) 



(46) 



From this one deduces that the optimal choices fox h — 2-k/p 
and Xm 



Mh are so that 



e 3 



I.e., 

h = 21/^1/3 {M + l)-2/3 ^ 1.58 (Af + l)-2/3. 
The corresponding estimated error becomes 



e{M) 



g-i2'»/3^(A/+l) 



-2.64 



(47) 



(48) 



I.e., if one wants to compute the normalization integral to 
P decimals precision one has to choose 



M+1 = ^P«0.87P. 

24/37r 



(49) 
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Fig. 3 This figure displays the computer time (measured in mil- 
liseconds) required to compute some normalization integrals to P 
decimals estimated precision. The cases considered are (i) the ground 
state wavefunction of the x*-potential {Eo), (ii) the 100* excited state 
wavefunction of the a;^-potential (i?ioo), and (iii) the ground state 



wavefunction of the (x 



-potential (Sww) with s — 1/100. 



C. Excited state wavefunctions of x'^ potential 

The main effect of having highly excited states lies 
in the extra factors in the error estimates, with C found 
in equation (39). Clearly C becomes large for large N . To 
account for this factor we must replace 10^^ by C^^ 10^^ 
in the error analysis. I.e., make the replacement 



P log 10 ^ F log 10 + TT tan 



N ■ 



(50) 



in the expressions like (49). There are also £'-dependent 
corrections to the remainding factors, but they are best treated 
numerically. 

D. Ground state wavefunction of {x'^ — 1)^ potential 
We next consider the lowest even parity eigenstate of 

-s^^A" + (x^ - 1)2^ = eV, (51) 

with s is small and positive. By WKB analysis one finds that 
the solution behaves like 



ip{x) 



-{x-iy{x+2)/3s 



(52) 



(up to an algebraic prefactor) in the region of interest (Rea; > 
0, |(a; — 1) (a; + 2)/3s| ^ 1). We again estimate Fourier 
transform of ■i/j{x)'^ by the saddle point method, assuming the 
relevant saddle point to be in a region where the approximation 
(52) is valid. The saddle point equation. 



(j)'{x,p) = 



d 

dx 



lf{x + 2)-ipx 



= 0, 



has a solution consistent with this assumption, 

-s = (l+if)^/^ 



The corresponding exponent is 

4 
3s 



(l + ips/2)3/2 - 1 



3s 



(53) 



(54) 



i)M'/^ 

(55) 



where the last equality is valid for ps ^ 1. To evaluate the 
integral to P decimals precision one must choose p = 2n/h 
so that (p{xs,p) = —P log 10 (or smaller). For very large P 
one may use the last approximation in equation (55) to find 

.,1/3 1.73 si/3p-2/3^ (56) 

(3 log 10)^/^ 

otherwise it is simplest to determine h numerically. Further, 
to evaluate the integral to P decimals precision one must add 
contributions from the (positive) x-range where 

« e-2(-i)^(-+2)/3s > 10-^. (57) 

I.e., we must take a^^ax to be the largest solution to the equation 
(x - \ f(x + 2) = (3/2) log 10 sP. For very large P the 
solution becomes 

Xn^ax « (flog 10)'/' S^I^P^I^ « 1.51 .1/3 (53) 

Hence the number of required function evaluations again grow 
asymptotically for large P like 

M+ 1 = w 0.87P, (59) 
h 
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cf. equation (49). 
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Fig. 4 These figures display the difference between the actually 
obtained precision Pobt and the estimated one Pest, for some 
normalization integrals, with precisions measured in decimal digits. 
The cases considered are (i) the ground state wavefunction of the 
a;*-potential (Eo), (ii) the lOO* excited state wavefunction of the 
a;*-potential (iJioo), and (iii) the ground state wavefunction of the 
(a;^ - l)^-potential (Bww) with s = 1/100. 



